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NOMENCLATURE 


a 0 (UJ) 

V'-J) 



0 

d 

dA 


local sonic speed 

sonic speed at freestream 

area of the freestream starting point grid 
element {x Q (I,J)} 

area of the surface impingement point grid 
element {"x m ( I , J) } 

least-square flow velocity interpolation coefficient 
array associated with the surface cell (i,j,k) 

buoyancy force acting on the droplet 

characteristic dimension of body 

droplet drag coefficient 

incompressible sphere drag coefficient 

fraction of total LWC contributed by Langmuir-0 
multi-droplet droplets of diameter d. 

fluid dynamic drag force acting on the droplet 

droplet diameter 

infinitesimal impingement area on body surface 

infinitesimal freestream cross-section corresponding to 
dA 


vi 


NOMENCLATURE 


d i ith Langmuir-0 multi-droplet droplet diameter 

[f-1 parametric cubic blending function basis; i=l,2,3,4 

G Cunningham drag correction factor 

G gravity force acting on the droplet 

g gravitational acceleration 

h stepsize used in the numerical integration of the 

trajectory differential equation 

i,j,k computational x,r,0 mesh node indices 

A A A 

i,j,k unit vectors // to x,y,z axis 

K(R v ) Stokes' parameter of droplet = c p( R v ) ‘ R v / 24 

K n Knudsen number = X/d 

T position vector along a trajectory segment 

L length of a trajectory segment = Ip - qj 

LWC liquid water content of droplet cloud at freestream 

LWC' liquid water content at body surface 

M Mach number of air flow restive to droplet 
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NOMENCLATURE 


CF 


MVD 


500 

p.q 

IQ I 

P. 


max 


r n * r t 


CF 

r(u,v) 


Mach number of air flow at compressor face 

Mach number of air flow at freestream 

mean volumetric diameter of the droplet cloud 

unit normal vector at body surface (pointing outward) 

unit vector // to V 

00 

p 

inertia parameter of droplet = p*V oo d /(18uC) 

freestream static pressure of air 

trajectory segment end points 

bi-cubic patch boundary condition matrix 

Reynolds number of air flow based on d 

Reynolds number of air flow relative to droplet 

radial boundary for cylindrical computational domain 

normal or tangential resiliency coefficient of 
particle 

inlet fan radius 

position on a bi-cubic parametric patch surface 


solid 
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NOMENCLATURE 


r u ( u » v ) 


r uv< u - v > 


Sao 


[U| 

[VI 

u 


U c’*c 


tangent vector along a constant-v parameter curve on a 
patch surface = 3r(u,v)/3u 

2 _ 

twist vector = 3 r(u,v)/3u3v 

total surface area of body surface (or a patch) 

arc length along a constant-9 curve on body surface 

freestream static temperature of air 

time, dimensionless with C/V 

00 

time 

2 3 

u-parameter basis function = [1 u u u ] 

2 3 

v-parameter basis function = [1 v v v | 
droplet velocity, dimensionless with 
droplet velocity 

potential flow velocity, dimensionless with V 
potential flow velocity 
freestream air velocity 

droplet velocity or Dosition computed at a corrector step 
(Adams predictor-corrector) 


lx 



NOMENCLATURE 



P 


w 

w ( i , j , k ) 


x 

x,r f 0 

x,y,z 

x Q (I.J) 

^( 1 , 0 ) 

x c (Uj) 

ct 

8 

Y 


droplet velocity or position computed at a predictor step 
(Adams predictor-corrector) 

normalized line (trajectory segment) length parameter 

sub-cell volume of a cylindrical field cell, with a 
corner at node (i,j,k) 

compressor face x-boundary 

freestream x-boundary of cylindrical computational domain 

droplet trajectory position, dimensionless with C 

cylindrical coordinates 

Cartesian coordinates 

freestream starting point grid 

impingement point grid 

impingement centroid grid 

angle of attack (pitch angle) 

local impingement efficiency 

ratio of specific heats of air = 1.4 

incident angle at impact (solid particle) 


x 



NOMENCLATURE 


Kronecker delta 

flowfield resolution (error) 

density ratio = p/p* 

integration error estimate at Adams predictor-corrector 
step 

maximum absolute discrepancy between U and U 

c p 

maximum relative discrepancy between U c and U p 

density of air 

density of droplet 

absolute air viscosity 

roll angle 

yaw angle 

levi-Civita anti-symmetric tensor 

angle between a point on sphere and xz-plane 

unit direction vector along the line from p to q 

freestream cross-section of the limiting envelope of 
trajectories 
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mean free path of air 
stepsize adjustment factors 
velocity relaxation time of a droplet 
Subscripts and superscripts 

derivative of a quantity 
quantity at freestream or origin 
nth iterate 
matrix transpose 
normal component 


tangential component 
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The extent and local flux rate of water droplet impingement on affected 
aircraft surfaces constitute the basic information needed for the design 
and analysis of various ice protection systems. 

The cause of aircraft ice accretion is mostly due to the presence of 
atmospheric clouds containing supercooled water droplets. The water 
droplet content of clouds generally decreases with altitude and beyond 
about 22,000 ft above sea level, called the icing altitude, clouds 
consist mainly of frozen particles and do not pose, in most 
circumstances, a significant icing hazard [1). For a given condition 
within the icing envelope smaller supercooled droplets may freeze 
entirely upon impact with the aircraft surface (rime icing), while 
larger droplets, requiring larger amounts of latent heat removal, may 
freeze slowly with runback (glaze icing). 

As the supercooled droplets impact on the surface, the governing 
transport parameter is the local droplet mass flux rate at the surface, 
which is in turn related to the normalized local surface flux function, 3 
(local impingement efficiency): 

Local droplet imoingement intensity = B-V. • LWC J ^f lref ^unit time 

Total droplet imoingement intensity = Voq-LWCJ" BdS 

S 

where LWC = liquid water content at freestream 

S = total surface area of body [unit area] 
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The definition of 8 as the local droplet flux rate normalized to the 
freestream flux rate follows from the continuity of droplet mass flow 
applied to an infinitesimal droplet stream tube (Figure 1) of 

differential area vectors dA at freestream (dA is not // to 

V in general) and dA at the surface of impact: 

(LWC) V w - dA^ = (LWC) V * dA (mass flow continuity) 

(LWC) n^- dA m = (LWC) ' ( -V*n) dA 

(LWC) mass flux at dA 

= (1.0) 

( LWC ) V mass fl ux at dA 

' CO OO 


n *dA 

oo c 


where 


n = unit normal vector at dA, 

n = unit vector // to V , 

00 ' 00 * 

V {\IJ = local droplet velocity vector at dA (dA^), 
(LWC) 1 = liquid water content at dA. 



Figure 1 - Oroplet Stream Tube 
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The ice protection systems currently in use or in development on various 
aircraft can typically be categorized into two kinds: 



System 

Source 

Typical Aircraft Application 

(a) 

anti -icing 

hot air (bleed) 

engine inlet, wing l.e.,ram 
air scoop, pitot tube 

(b) 

de-icing 

pneumatic boot 

wing 1 .e. , rotor blade 



(pulsed air) 
electrothermal [ 2 1 

wing l.e., pitot tube, rotor 
blade, stabilizer 



electroimpulse [3] 

wing l.e., rotor blade, 
strut, stabilizer 



fluid (freezing 
point depressant) 

wing l.e. [4| 


In an anti-icing system, heat is continuously supplied (during system 
operation) to the affected surface, so that all of the impinging water 
can be evaporated or maintained above freezing. Thus, for the steady- 
state heat transfer analysis of the system, the total as well as the 
local water impingement intensities must be known at the icing 
interface. 

As for the de-icing system, ice is allowed to accumulate to a certain 
level and heat or some form of mechanical energy is supplied in a 
transient manner to shed off the ice, thereby saving substantial energy 
expenditure (compared to the anti-icing system) at the expense of the 
aerodynamic penalty due to ice accretion. For the performance analysis 
of a de-icing system, the extent and shape of the ice accreted must be 
known as well as the tolerable level of ice accretion in terms of the 
associated aerodynamic penalty. Analytically, this requires ice 
accretion modeling [5,61 » for wnich the local water impingement 
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intensity is an input, and detailed flow analysis about the body with 
the prescribed ice shape [7,81. 

Thus the determination of local water droplet impingement efficiency (8) 
on aircraft surfaces is a fundamental task in quantifying the aircraft 
icing phenomena. 

Oespite its importance in the design/analysis of ice protection systems 
for most of the present day engine inlets and wings which are highly 3- 
dimensional, there has been very little analytical or experimental work 
to determine the water droplet impingement efficiencies on 3-dimensional 
configurations. 

The purpose of this research work is to develop a 3-dimensional particle 
trajectory analysis computer code to predict the local water droplet 
impingement efficiency (B) on a representati ve commercial turbofan type 
engine inlet. This work grew out of the need to develop analysis tools 
leading to improved engine anti -icing and sand separator systems at the 
Boeing Military Airplane Company (BMAC), Wichita. 
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2.0 PREVIOUS RELATED RESEARCH 


A detailed review of the relevant research literature is discussed in a 
recent report by Shaw [9|. 

The bulk of the available literature in the area of B determination 
comes from the extensive NACA icing research efforts in the 1940-1950 
time period. The NACA research program concentrated mainly on the 
experimental determination of B on axisymmetric geometries [10,11] and 
2-D airfoils ( 12-16]. 

One of the earlier analysis efforts is due to Langmuir and Blodgett [17] 
who calculated water droplet trajectories to predict impingement 
efficiencies about circular cylinders using a differential analyzer. 

More recently, a number of researchers developed several water droplet 
trajectory codes to compute 8 [18-20] as well as to model ice accretion 
[5,6] and to assess aerodynamic penalties [7,8] on 2-D airfoils. Code 
development applicable for engine inlets was limited, partly due to the 
more complex flowfields involved and due to a complete lack of test data 
on these geometries. However, a code [21,22] was developed by BMAC and 
applied to various axisymmetric engine inlet anti-icing analysis 
problems. 

Trajectory code development for 3-dimensional impingement problems has 
been of limited extent. Although two codes [23,24] exist that are 
capable of analyzing 3-D impingement problems, calculation of B was not 
reported . 

The 3-0 code developed by Norment [23] uses the Hess-Smith 
incompressible panel potential code [25] and a variable order Adams 
predictor-corrector integrator to solve the trajectory differential 
equation. To compute the fluid dynamic forces acting on the droplet at 
each trajectory position, it uses the direct approach of computing the 
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flow velocity at the required position by summing over all the panel 
source and vorticity contributions. Since each droplet trajectory 
involves computing hundreds of intermediate trajectory steps, the 
computing time will generally be high. As with other panel flow codes, 
the accuracy of computed flow velocities is ultimately limited by the 
panel density. For a trajectory segment near the panel surface, the 
intermediate trajectory steps will be crowded (compared to a segment far 
away from the surface) because of stronger flow gradients there. If the 
mean distance between these intermediate steps is small compared to the 
linear dimension of surface panels, the direct approach will result in 
non-smooth flow velocities along the trajectory and lead to numerical 
problems in solving the trajectory equation of motion. Thus, the choice 
of the direct approach for the trajectory code seems questionable in 
view of the inherent danger and high computing times involved. This 
code currently uses the non-lifting version of the Hess-Smith panel 
code, and is not applicable for problems involving engine inlets. 

The recent 3-D code developed by Stock [24) employs a finite volume 3-D 
Euler flow code [261 and a 4th order Runge-Kutta scheme to solve the 
trajectory equation of motion. It was applied to the droplet 
impingement problem on a 3-0 engine intake, utilizing a body-fitted 
computational mesh (grid approach). Body-fitted mesh definition of 
computational flow domain is generally accepted as one of higher flow 
resolution than any other fixed orthogonal mesh systems, because of its 
grid adaptability near the boundaries [ 27 | . However, because of the 
finite volume approach employed, uniform flow velocity was assumed 
throughout the volume of each mesh cell, while at least five 
intermediate trajectory steps were computed in each mesh cell. This 
approach may be acceptable in the far field region, where the flowfield 
is approximately uniform. Near the boundary surface, the assumption of 
flow uniformity can be incorrect in predicting particle impact on the 
surface. The tangent impact points computed from this code indeed 
reveal erratic jumps between several pairs of adjacent tangent 
trajectories [24|. 
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Since it is the fluid dynamic forces acting on the droplet that 
determine the droplet trajectory, accurate flowfield definition is a 
prerequisite for accurate droplet trajectory computation. For the grid 
approach, the additional requirement of accurate flow velocity 
interpolation along the trajectory must be met. Also, an accurate 
surface geometry definition of the body is needed in order to locate the 
impingement points precisely. This is not a trivial task for 3- 
dimensional geometries such as engine inlets. B is a local surface flux 
function and its accuracy is very sensitive to the local surface 
geometry. Flat panel surface definition, as is done in panel potential 
flow analyses, will not be adequate unless sufficient panel density is 
used. 

For the present investigation, the grid approach is adapted based on the 
consideration of computational efficiency. Flow velocities were 
computed using a 3-0 compressible full potential flow code [28,29) on a 
cylindrical mesh system. Linear and least-square interpolation 
techniques are employed for flow interpolation along trajectories to 
ensure smooth and accurate resolution of the flowfield. State-of-the- 
art bi-cubic parametric surface modeling techniques [30| are utilized to 
obtain an analytical definition of the 3-0 engine inlet surface studied 
as well as to compute the impingement points accurately. Variable step 
fourth order Runge-Kutta and Adams predictor-corrector integration 
schemes were used to solve the trajectory equation, together with an 
automatic stepsize control scheme to maintain the desired integration 
accuracy in the numerical solution of the trajectory differential 
equation. 
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3.0 TRAJECTORY MODEL 


The general motion of droplets moving through turbulent air flow regimes 
is not considered in this study. A rather simplistic approach, taken by 
researchers as early as 1940' s [17], is to describe the quasi-steady 
motion of small spherical droplets moving in the steady flow of air, 
while the motion of droplets does not disturb the air flow. The 
predominant force acting on a droplet is then the fluid dynamic drag 
arising from the relative (slip) velocity of air with respect to the 
droplet. This is a valid approach in view of the fact that for the 
typical icing design conditions of intermittent and continuous maximum 
[ 1 1 , the maximum concentration and mean volumetric diameter (MVD) of 
droplets are: 

intermittent maximum continuous maximum 

LWC MAX a 3.0 gm/m 3 LWC MAX * 0.8 gm/m 3 

MVD MAX = 50 pm MVD MAX * 40 pm 

For the concentrations and sizes of droplets within the icing envelope, 
the assumption of undisturbed airflow and spherical shape (due to 
surface tension) of droplets are quite valid. 

3.1 Model Assumptions 

(1) Single phase (air) flow about the body - particle phase does not 
disturb the flowfield of the gas phase. 

(2) Quasi-steady-state approximation - at each instant and position, 
the steady-state drag and other forces act on the particle. 

(3) Compressible or incompressible potential flowfield of the gas phase 
about the body. 



(4) Sphenca] shape of particles. 


Additionally, viscous flow effects such as thick boundary layer 
formation and flow separation are not considered because particle 
impingement usually occurs in the forward part of the body. 

3.2 Trajectory Differential Equation 

Under the model assumptions, the forces acting on the particle are the 
fluid drag, buoyancy, and gravity. By applying Newton's second law and 
non-dimensional i zing (Appendix B), the particle equations of motion 
reduce to the following: 

dU i /dt = C 0 (R v )-R-(V r U.)/(24P) - (l-o)g CS^/V* (i = 1,2,3) (3-1) 


where 

* 2 

P = p d /(18uC) = inertia parameter of droplet, 

t £ time (dimensionless with C/V ), 

GO 

a e p/p* = density ratio of air to particle, 

C e character ist ic dimension of body, 

R v e relative Reynolds number of droDlet, 

U e particle velocity (dimensionless witn Vj, 

V e potential flow velocity (dimensioness with V ). 

Because of the way the slip velocity, V - U, appears in the slip 
Reynolds number (R ), equation (3-1) must be solved numerically in 
general: in some ideal cases, when V is a simple function of position 

and R can be expressed in a special form, equation (3-1) can be solved 
analytically (31, 32 I - 
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3.2.1 Meaning of P (inertia parameter) 


For a particle injected into the uniform flow, V = V , and obeying the 
Stokes' law of drag (Cg = 24/R^): 

dU/dt = (V^-UJ/P (3-2) 

(neglecting (l-oJgC/V^ << 1.) 

Equation (3-2) can be written as 

dOJ-VJ/dt = -(U-V^J/P (since is constant) 

which integrates to 


U = V o - (U Q - VJ exp(-t/P) ; (U Q = U(t=o)) (3-3) 

From equation (3-3), U ■+ 7^ as t •» » monotonical ly , i.e., velocity of 

the particle relaxes to the flow velocity after a long time. Thus P is 
the non-dimensional equivalent of the velocity relaxation time (t ) 
characteristic of the particle: 

P • \(v c ) 

x v = p d^/(18p) [unit timel (3-4) 

Equation (3-4) implies that the larger, heavier particles will take 
longer time to relax to the flow velocity than the smaller, lighter 
particles. In the general case of arbitrary flow and Cg t Cg (Stokes), 
the velocity relaxation time concept is still useful in that a rough 
order of magnitude estimate of the particle motion in a given flow can 
be obtained. 
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3.2.2 Scaling of the Trajectory Problem 


As was done in the previous section, the trajectory equation (3-1) can 
be simplified for the case of Stokes' drag law: 

dU./dt = (V i -U.)/P - (l-o)qC6. z /\ll (3-5) 

Equation (3-5) implies that, neglecting terms due to gravity and 
buoyancy, the particle trajectory for a specified starting condition is 
completely determined by the inertia parameter P for all dynamically 
similar flows. The only trajectory similarity parameter required is P, 
in addition to the usual flow similarity parameter such as the Reynolds 
and/or Mach number (depending on the degree of flow compressibility). As 
long as the Stokes' law of drag holds along the trajectory, matching of 
P guarantees the trajectory similarity for the same set of initial 
conditions for all dynamically similar flows. This similarity concept 
for constant P breaks down at the limit of the Stokes' law of draq 
because K(R ) is non-linear in R y ; its deviation from unity at a point 
in the trajectory is a measure of the extent to which the drag 
coefficient differs from the Stokes' law value. Consequently, the 
trajectory problem cannot be scaled in general due to the trajectory 
dependent Stokes' parameter, K(R v ). 

3.3 Drag Coefficient for Spherical Particles (Cg) 

The particular form of the drag coefficient used in this study 
incorporates an analytical form for the standard drag curve and the 
Cunningham drag correction for molecular slip and compressibility 
effects: 

C 0 (M,R V ) = Cg lnC ’(R v )/G(M/R v ) (3-6) 


where 


Cg inc ‘ = incompressible sphere drag coefficient 
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Cunningham drag correction factor 


G(M/R y ) = 

3.3.1 Incompressible Sphere Orag Coefficient (Cg 1nc ‘) 

There exist some experimental drag data [33,34] on water droplets in 
sizes well above the millimeter range where the droplets tend to deviate 
from sphericity. The effects of droplet instability/break-up are 
pertinent for larger droplets and hence are not considered here. 

The deviation of water droplet drag data from that of the standard drag 
was found to be significant only for droplets of diameters larger than 
about 1 mm and for Reynolds numbers greater than about 1000. This 
observation was mainly due to the flattened shape of droplets in the 
size and Reynolds number range studied in the still air settling speed 
measurement [ 33 ] . In fact, a recent investigation [35] reported that, 
in both the Navier-Stokes flow analysis and settling speed measurement 
results, no significant differences in drag larger than the measurement 
errors were found between the solid and liquid spheres. 

Equation used is the integrable form of Putnam [36], 

C D inC ’(R v ) = C 0 Stokes (R v ) • (1 + g R v 2/3 ) (3-7) 

which agrees to within about 5 % of the standard drag curve in the range, 
0 < R v <1000. Comparison of this equation with several other available 
forms is listed in Appendix C; equation (3-7) is listed as C^ (Putnam). 

3.3.2 Cunningham Orag Correction (G(M/R v )) 

For small droplets less than about 5 urn diameter, reduction in drag can 
occur because of the molecular slip of air. Whenever the size of the 
particle becomes comparable to the mean free path of air molecules, this 
non-continuum effect can be significant. The first attempt to correct 
for this was made by Millikan in his oil drop experiment. He used the 
following correction formula to the Stokes' viscous drag for oil 
droplets: 
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C 0 = C 0 Stoke 5 (1 + (X/d) • (C 1 + C 2 e C3dA )j- 1 (3-8) 

where X is the mean free path of air, 
d is the particle diameter, 

C^, C 2 and are empirical constants. 

The factor X/d is also known as the Knudsen number (K ) which can be 
shown to be proportional to M/R v from the kinetic theory of gases [37): 

X/d = K n cc M/R v 

where M = Mach no. of gas flow relative to the particle. 

The form of the correction adapted in this study is due to Cal son and 
Hoglund [381, who proposed the following empirical fit to available 
experimental data for the ranges M<2.0 and R v <1000: 

G(M/R y ) - A/B (3-9) 


where 


A i 1 + (M/R y ) [3.82 + 1.28 exp (-1.25R y /M) 1 , (3-10) 

8 = 1 + exp (-.427M" 4 ’ 63 -3R y "- 88 ). (3-11) 

The numerator in equation (3-9), A, has the same form used by Millikan 
and only the numerical constants have been modified. This term 
reoresents the drag reduction factor to the incompressible drag due to 
the molecular slip or rarefaction effects. 

The denominator, B, in equation (3-9) is the additional correction to 
account for the Mach number dependence of the particle drag 
(compressibility) in continuum flow. 
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It must be noted that the Cunningham correction, expressed as in 
equations (3-9), (3-10) and (3-11), must be evaluated at every 

trajectory step in the flowfield. The relative Mach number, M, can be 
evaluated from the compressible Bernoulli equation applied to the 
potential flow velocity: 

vf /2 + a*/( Y -l) = V 2 /2 + a 2 /(y-l) (3-12) 

where a 5 /yRT $ , 

a = / W, • 

00 1 5 00 

Thus, 

T s ■ 2VR (1 - O'/VjVv** (3-13) 

M = |V-U| • MJSy RT s (3-14) 

Substitution of (3-13) into (3-14) gives 


m = iv-ui • [ (y- 1) (i-v 2 )/2 + y rt sco • v; 2 r* 


(3-15) 
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4.0 COMPUTATIONAL METHOD 


Analytical determination of water droplet impingement efficiency 
involves calculation of the flowfield about the body and calculation of 
individual particle trajectories that lead to impingment on the body. 
The local impingement efficiency (8) is obtained by computing the 
impinging surface droplet flux relative to freestream flux as a function 
of body surface coordinates. 

4.1 COMPUTATIONAL PROCEDURE 

Droplet impingement analysis on the 3-D engine inlet involved the 
following major steps: 

(1) Bi-cubic parametric description of the inlet surface. 

(2) Potential flow analysis about the inlet. 

(3) Numerical integration of the trajectory equation. 

(4) Calculation of the limiting envelope of trajectories. 

(5) Calculation of 8 from the intermediate trajectories within the 
limiting envelope. 

The steps involved in the procedure are illustrated in Figure 2. 

4.2 8 1 -CUB I C SURFACE PARAMETRIZATION 

Any point on a 3-dimensional surface element (patch) can be analytically 
defined in terms of a set of patch corner boundary conditions through 
bi-cubic surface parametri zation (See Appendix A). Parametrization is 
complete when all the patch corner boundary conditions (patch boundary 
matrices) are obtained for the particular system of patches making up 
the composite patch surface of the inlet. This procedure involves cubic 
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FIGURE 2 - COMPUTATIONAL PROCEDURE 



parametric splining of the two sets of curves defining the composite 
path surface, using the accumulated chord length parametrization and 
Hermite interpolation schemes [39 j, to compute the required end point 
tangent and twist vectors for each curve segment. 

The purpose of parametric surface description is two-fold: 

(1) It provides accurate surface normal velocity boundary conditions 
required for the full potential flow code input. 

(2) Accurate trajectory-surface intersections (impingement points) can 
be obtained through such parametrization (Appendix A). 

A wiremesh diagram of the inlet derived from the bi-cubic patch 
parametrization is shown later in Figure 15. In the figure, the 
straight line edges of the wiremesh patches are for illustration only. 
These do not represent the actual patches whose edges are curved. For 
the inlet investigated in this study (737-300 prototype inlet), about 
600 patches were used to define the surface. 

4.3 POTENTIAL FLOW ANALYSIS 

Flow velocities are computed by the 3-D full potential code [29 1 on 
cylindrical mesh grids (69 x meshes, 49 r meshes, and 16 9 meshes). The 
flow code solves the full partial differential equations of compressible 
transonic potential flow by a finitq difference scheme. The convergence 
acceleration is achieved by the successive line over-relaxation (SLOR) 
and multigrid techniques [ 40 , 4 1 J . The multigrid scheme utilizes four 
levels of coarse and fine grids about the original mesh chosen such that 
during iteration cycles flow solutions are passed from one level to 
another to achieve the extremely fast convergence of flow solutions. 
For an average engine inlet flow problem involving 50,000 mesh points, 
the CRAY-IS computing time is only about one minute. 

A typical adaptive mesh grid used for engine inlet flow analyses is 
shown later in Figure D8, Appendix D. 
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4.3.1 Flow Accuracy (3-0 full potential code) 


This particular CFO (Computational Fluid Dynamics) computer code used 
for this investigation is an example of a time-tested code. Since its 
initial production version developed in the early 1970s, many validation 
comparisons with experimental data are available, including the NASA 
wing-pylon-nacelle model (Figure Dl, Appendix 0). Practically all of 
the Boeing commercial engine inlets were designed using this code and 
later correlated with wind tunnel data. However, most of these 
validations are Boeing proprietary data and cannot be included here. 
Figures 02 through D7 (NASA wing-pylon-nacelle model) and Figures 09 
through Dll (full scale commercial turbofan-engine type inlet) are from 
NASA CR3514 and these show good agreement with the measured surface Mach 
number data. 

4.3.2 Flow Velocity Interpolation at Trajectory Steps 
Two different interpolation schemes are employed: 

(1) Volume weighted linear interpolation in field cells. 

(2) Least-square interpolation in surface cells. 

Example of the two types of mesh cells are illustrated in Figure 3 
showing a typical coarse mesh definition of the flow domain about an 
engine inlet. 

4.3.2. 1 Linear Interpolation Formula 

Figure 4 depicts a cylindrical field cell with the mesh node origin at 
(i,j,k); i, j and k are the x, r and 9 mesh indices of the node. The 
flow velocity, V(x,r,9), at an interior point, (x,r,9), is 
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V(x,r,e) = (Vol)" 1 - {V(i,j,k)-w(i+l,j+l,k+l) + V(i , j,k+l) -w(i+l, j+l,k) 

+ V(i+l,j,k+l)-w(i,j+l,k) + V(i+l,j,k)-w(i, j+l,k+l) (4-1) 

+ V(i , j+l,k) -w(i+l, j,k+l) + V(i , j+l,k+l) -w( i+1, j,k) 

+ V(i+l,j+l,k+l)-w(i,j,k) + V(1+l, j+l f k) -w(i , j,k+l) } 

where 


‘i i x i x i + i- r j i r i r j + l- k+ l- 


V(i,j,k)= Flow velocity at the node (i,j,k). 


vol E ( x ,+r x i > • ( r j+r r j ) ' (e k+r 9 k > /2 

= total volume of the cell, 

w(i,j,k)= volume of the sub-cell whose two corners are 
(x,r,e) and (x i ,r j ,0 k ). 


Explicitly, w's are: 

w(i,j,k) = (x-x.)-(r 2 -r 2 )-(0-0 k )/2 

w(i+l,j,k) = (x - +1 ~x) • (r 2 -r 2 ) • (0-0 k )/2 

w(i,j+l,k) = (x-x i )-(r 2 +1 -r 2 )-(0-0 ]< )/2 

w(i,j,k+l) = (x-x i )*(r 2 -r 2 )-(0 k+l -0)/2 

w(i,j+l,k+l) = (x-x i )-(r 2 +1 -r 2 )*(0 k+1 -0)/2 

w(i+l,j+l,k) = (x i+1 -x)*(r 2 +1 -r 2 )*(0-0 k )/2 
w(i+l,j,k+l) = (x. +1 -x)*(r 2 -r 2 )-(0 k+1 -0)/2 
w(i+l f j+l,k+l) = (x i+1 - x )-( r j +1 - r2 )*( 9 k+L -9)/2 
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i+1 ,j+l ,k+l ) 


FIGURE 4 - Cylindrical Field Mesh Cell 
4. 3. 2. 2 Least-Square Interpolation Formula 

When a mesh cell intersects the surface, as in the surface cell shown in 
Figure 5, interpolation becomes difficult. For example, a straight 
forward application of the Taylor series formula would require 
evaluation of the flow gradients with respect to the coordinate 
variables which depends on the particular way the surface intersects the 
mesh cell. This is a time consuming process since all the possible 
cases of surface intersecting cell geometries have to be accounted for. 

A different approach is employed in this study, whereby the flow 
velocity at a desired point is assumed to be a function of the space 
coordinates and the unknown set of coefficients are to be determined 
from the least-square fitting of this function at the exterior and mesh- 
surface intersection points (Figure 5) associated with the surface cell: 
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(4-2) 


V(x 1 ,x 2 ,x 3 ) = a 1 + a ? x 1 + a 3 x ? + a4 x 3 + a 5 x l x 2 + a 6 x 2 x 3 + a 7 x l x 3 

This is equivalent to the lowest second order Taylor expansion about the 
point (x°,x°,x°) in a vector form: 

V(x r x 2 ,x 3 ) = 7(0) + Iri^Vj^k^x^o + 0(3) 

= 7(0) + + x 2^s7 2 ^o + x 3^ 3 ^o 

+ l/2[ x 1 x 2 (-|^L-) 0 + X 2 X 3^ 3x~3x 3 ^o + X 1 X 3^ 3x^3x~^ 


+ 0(2;higher order) + 0(3) 


(4-3) 


where (0) and subscript o are taken to mean the quantities evaluated at 
some arbitrary origin (x°,x 2 ,x°) in the neighborhood of the surface 
cell. 


Thus the least-square model, equation (4-2), reduces to the lowest 
second order Taylor expansion of the flow velocity if the following 
equalities are made: 


a 1 = V(0) 


*i ■ 1 ■ z - 3 ' 4 


a 5 - 1/2 ( 


a 6 = 1/2 ( 


a ? = 1/2 ( 


2 - 
3 c \l 


3x^3x 2 'o 


2— 

3 


3x 2 3x 3 'o 


2— 
3 v 


3x^3x 3 'o 
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Equation (4-2) can be put in the matrix form 


U, 5 


1 s 2 *3 *4 *5 *6 *7 ^ 


7J 


( V ( x ^ , X£ , x^ ) ] 


(4-4) 


where 


= 1, ^2 ~ x ^ > £3 = x 2 * ^4 x 3’ ~ x i x 2» £g = x 2 x 3 ~ x ^ x 3' 


For a total of n (n>7) boundary points (exterior mesh and surface-mesh 
intersection points) defining the surface cell, the successive 
application of equation (4-4) to the n boundary points results in: 


*1,1 *1,2 *1,6 *1,7 




r— - 

V(l) 

*2,1 *2,2 *2,6 *2,7 

*3,1 

• • • • • 


*1 

*2 

*3 


7(2) 

• 

• 

• • • • • 

« • • • • 

• • • « • 

* • • • • 


d 4 

*5 

*6 

. * 7 


V(n-l) 

|_ s n, 1 *n,6 s n,7_j 




7(n) 


This is an over-determined system with 7 unknowns (a's) and n equations. 
V(k) is the known flow velocity at the k L boundary point. The matrix 
[£(I,0)| ( 1 = 1,2, ... ,n; J=l,2,..,7) contains the terms involving only the 
coordinate positions of the n boundary points defining the surface cell. 
Equation (4-5) is solved for a's using Householder's least-square 
minimization procedure for over-determined system of equations ( 42 | : 
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E = IU C lm a m - V 1 )‘ 

l=lm=l im m 1 


; 1 = 1,2, ,n 

m = 1 ,2, ,7 


Minimizing E with respect to ( k = 1,2,' 

°-* ■ ,LU [(5 ^- r ' )2] 

"l=L=l ‘ V-l 

= < hA - “i^rn/W 

- 2 j/«lk) T -<«lk r k * f l> 


Therefore [sfCOOT = [s] T [v] 


where [£] = [£ lm ]. 


[A] = [a 1 ^2 a^ a^ a^ ag a^] » 


[V] 5 [v 1 v 2 


V , V ] T . 
n-1 n J 
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In some special situations when n < 7, extra mesh or surface-mesh points 
(and associated flow velocities) adjacent to the surface cell are u^ed 
to meet the matrix dimensional requirement, n > 7, in solving for the 
least-square coefficients (a^. }. 

To achieve computational efficiency, all of the surface intersecting 

mesh cells were identified and the least-square coefficients determined 

and stored in a file prior to the running of the 3-D trajectory code. 

Each mesh cell is identified, as shown in Figure 5, with the set of mesh 

node indices, (i,j,k), and the associated set of 21 least-square 

coefficients, (a^(i,j,k)}, 1=1, 2, 3; m=l,2,...,7 (1 is the coordinate 

component index for computing V.). As long as a trajectory position 

* 1 

(x ,r , 9) is inside the surface cell (i,j,k), the coefficients (a (i,j,k)} 
are used in equation (4-2) to compute the interpolated flow velocity at 
that position (x=x^, r=x,,, 0=x^). 

A= mesh surface 
intersection 
point 

0= interior mesh 
node 

• = exterior mesh 
node 

FIGURE 5 - 2-0imensional Illustration of a Surface 
Cell Associated with an Exterior Mesh Node (i,j,k) 

For the surface cell (Figure 5) associated with the node (i,j,k), the 
least square coefficients (a (i,j,k)} can be obtained from the positions 
and flow velocities at 8, C, and the node (i,j,k). 

The least-square interpolation formula, equation (4-2), resolves the 
potential flow velocity near the surface very accurately and smoothly. 
This is shown in Figure 6 where the comparison is made between the 
interpolated and the CFO output velocities at the mesh-surface 
intersections on the lower cowl surface of the 737.300 prototype engine 
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RADIAL 


Figure 6 - Lower Cowl Surface Fl< 
Least-square Interpol. 





inlet. It should be noted that the interpolated curves are obtained not 
by fitting the CFD output flow velocities at the surface (solid symbols 
in Figure 6) but from the least-square coefficients of the surface cells 
involved, by using equation (4-2) continuously along the cowl surface. 

4.4 Numerical Integration of the Trajectory Differential Equation 


Substituting the drag coefficient, expressed as in equations (3-6), 
(3-7), (3-9), (3-10) and (3-11), into the trajectory equation (3-1) 


dU,/dt - C ^ tokes (R v )-R v -(l ♦ iR v 2/3 )-(V j -U f )-[24P-G(M/R v )r 1 (4-6) 

- a-d)s cs i 2 v 2 


where G(M/R v ) 


1 + (M/R^) ( 3.82+1.28exp(-1.25R v /M) ] 
1 + exp[-.427M -4,63 - 3R“ * 88 ] 


M = M(V,U) as shown in equation (3-15) 

R v = R v (V,U) as defined below equation (3-1) 
V = V(x) 


From the above functional relationships, the R.H.S. of equation (4-6) 
depends only on the particle position and velocity: 

dU/dt = F(x,U) 

Together with the definition of U as the time derivative of 7, we arrive 
at the following two coupled first order differential equations in time: 


dx/dt 

= u. 

(4-7) 

dU/dt 

= F(x,U). 

(d-8) 
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The above represents an initial value problem and thus requires a self- 
starting type procedure in the numerical integration. The Runge-Kutta 
method is a self-starting type with high degree of accuracy. One 

disadvantage is that a large number of function evaluations is involved. 

Compared to the same order Runge-Kutta integrator, the Adams predictor- 
corrector requires one-half the number of function evaluations with 
comparable accuracy, but is not self-starting. 

In this study, a combination of the 4th order Runge-Kutta and Adams 
predictor-corrector schemes with an automatic stepsize control is used 
to solve the trajectory equation (4-6). A flow chart of the numerical 
integration scheme is snown in Figure 7. 

4.4.1 Runge-Kutta Scheme (4th order) 

The 4th order scheme, accurate to 5th order in Taylor expansion, is used 
to start or restart the integration process from the initial condition 
or when the stepsize is changed due to the error control process at the 
end of an Adams predictor-corrector step. 

The coupled first order equations (4-7) and (4-8) take the following 
Runge-Kutta forms 

U(n+1) = U(n) + h* (A + 2B + 2C + D|/6 (4-9) 

x(n+l) = x(n) + h * { a + 2b + 2c + d(/6 (4-10) 

where 

h = Vl - V U(n) E U{t n } ’ *(") 5 7(t n } 

a = U(n), A = F[x(t n ),U(t n )| 

b = a + h-A/2, 8 = F[x(t n )+h-a/2,b| 
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c = a + h*(A+B)/4, C = F(x(t n )+h* (A+b)/4, cj 

d = c + h-C/2, 0 = F[7(t n )+h-(A+b+2c)/4, d| 

4.4.2 Adams predictor-corrector scheme (4th order) 


This scheme predicts and corrects the next time step (t n+4 ) from the 
three previous Runge-Kutta time steps (t n+ p fc n+2’ ^+3^* 


Predicted 


Corrected 


U p (n+4) = U(n+3) + h*[-9F(n) + 37F(n+l)-59F(n+2)+55F(n+3) ]/24 
x p (n+4) = x(n+3) + h-[-9U(n) + 37U(n+l)-59U(n+2)+55U(n+3) |/24 
U c (n+4) = U(n+3) + h- (F(n+l)-5F(n+2)+19F(n+3)+9F{x p (n+4) ,U p (n+4) } |/24 
x c (n+4) = x(n+3) + h- (U(n+l)-5U(n+2)+19U(n+3)+9U p (n+4) 1 /24 


where 


F(n) s F[ x(n) ,U(n) | . 

This procedure is recursive, i.e., as long as the agreement between the 
predicted and corrected values are within the specified error margin 
integration proceeds in a step by step manner. 

4.4.3 Automatic Stepsize Control 

After each predictor-corrector computation, the integration error (e) is 
checked to determine whether to accept the corrected values and proceed 
or to reject the step and restart the integration using the Runge-Kutta 
procedure. The integration error at a particular step is not the 
truncation error occurred at that step corresponding to the particular 
choice of numerical scheme, but is the global error of the numerical 
solution from the true solution at that step. The approach of 
controlling the truncation error at every step usually involves 
additional number of function evaluations comparable to those required 
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FIGURE 7 - INTEGRATION SCHEME FOR TRAJECTORY COMPUTATION 
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to integrate the equation, and still does not guarantee satisfactory 
control of the global numerical error. 


The approach adapted here is to use the discrepancy between the 

predicted and corrected values as a measure of the integration error, 

e = e(U c , U p ). The two types of error indicating the discrepancy 

between U and U are the absolute, e, (U^-U^ ) , and the relative 
pc 1 . c p 

e 2 ( U p/U<l), errors: 

Type of 

Situation Error Preferred Min. Error ( smin) Max. Error ( anax) 


(1) I U c I , I Up I large relative (e 2 ) 


0 . 01 ( 1 %) 


(2) I U c | , | Up I small absolute (z { ) 


SV(f1owfie1d error) 


The error bound on must be tied to 5 V, since the R.H.S. of the 
differential equation (4-6) depends on (V-U); if |V| is accurate tc SV, 
then the absolute error of U should be of the same order of magnitude. 
However, if is within the flowfield accuracy, e 2 may be unacceptably 
large; 


I I* I Up I ma y be sma11 (comparable to SV) but 

l(uj-uj)/ujl or I l-iuj/ujll may be large. 

Thus, the flowfield accuracy (SV) plays an important role in the error 
control process. SV for the 3-0 inlet flowfield was about 0.001 based 
on the flow velocities normalized to the Mach number at freestream 
(M„ = .267). 

Based on the above physical considerations, the following integration 
error estimate/control scheme is devised: 
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(1) 


e l’ £ 2’ and initialized to zero (i=l,2,3) 

(2) lujl < 6V ; ei (i) = luj - ujl 

IU*I > 6V ; e 2 (i) = MAX { I (uj-Uj )/u’ I , ll-luj/U*| I } 

(1=1, 2, 3) 

= MAX (e lf MAX le 1 (l) f e 1 (2),e 1 (3) |) 
e 2 = MAX (e 2 , MAX [e 2 ( 1) ,e 2 (2) ,e 2 (3) I ) 

(3) case 1 (e^ < 2 SV) ; 

emin < e = e^ < emax -*■ successful predictor-corrector step 

e = e 2 < emi n rejected, h' = h-X^ (increase stepsize) 

e = e 2 > emax -*■ rejected, h' = h/ X 2 (decrease stepsize) 

case 2 (e^ > 2 SV) -*• rejected, h' = h/X 2 (decrease stepsize) 


(4-11) 

(4-12) 


Thus, the error control scheme first checks to see if the maximum 
absolute error (e^) of (U c , U ) components is within the flow 

resolution. [f it is (case 1), then it checks whether the maximum 
relative error (e 2 ) is within the set relative error margin 

(emin, emax). Otherwise (case 2), the step is rejected, and stepsize 
decreased to restart the integration process using the Runge-Kutta 
procedure. 
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For the 3-0 inlet trajectory analysis, the following error control 
parameters are used: 


<5V = 

0.001 


emin = 

0.001 

, emax = 0.01 

x i - 

1.5 

, X 2 = 1.87 


Integral or integer fractional relationship between X^ and X^ is to be 
avoided because of the danger of repetitive stepsize changes. 

4.5 Calculation of the Limiting Envelope of Trajectories 


The limiting envelope of droplet trajectories is the surface traced out 

by the inner and outer tangent trajectories (Figure 8). For droplets 

starting at freestream constant-x plane with the same initial velocity 

as V , the bounding radial starting positions, (r min(9 ), r max(9 )), 
°o 0 0 0 0'' 

are searched for each selected value of 9 q which results in a pair of 
tangent trajectories. Repeating this process by sweeping 9 q with 

selected increments, the freestream impingement bound, r(r Q min(9 o ), 
r Q max(9 0 )), is determined. T represents the cross-section of the 
impinging envelope of trajectories at the freestream constant-x plane. 


A tangent trajectory is determined via a trial and error process. Along 

a radial line (9 Q ray) on the constant-x freestream plane, a pair of 

radial starting positions r Q ^ and r 2 are searched that result in 

impinging and non-impinging trajectories , respectively, to the engine 

cowl surface. Once r . and r _ are found, an iterative bi-section 

ol o i 

procedure is applied until a set tolerance is met: 


( 1 ) 

( 2 ) 

( 3 ) 


r Ql (impinging), r Q2 (non-impinging); (x q , 9 0 ,U 0 = VJ fixed 


r o = (r oi + r o2 )/2 (new 9 uess) 



ol 


I < TQl 7 


Y E S _ {tangent trajectory found} » ( £XIT) 
N 0 — ►( 4 ) “ ~ 
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( 2 ) 


( 4 ) 


compute trajectory 
until particle f 


impinges: set r Q ^ = r Q — 

re ach e s~compu t a t i ona 1 
boundary without impinging 


<5) 


(5) r o = r o2 ^ 

If r Q 2 had the trajectory end point radial position (at the fan face, 
the downstream computational boundary) within the fan radius, 

r min = r ,. Otherwise, r max = r 
o ol o ol 

Steps (1) thru (5) are performed for 0 Q (n) = 0 Q (n-l) + A0, 

n = 1,2, • • • ,nmax; nmax = 30, 0 0 (0) = 0° and A0 = 180°/nmax. 

The 2(nmax+l) tangent trajectories found this way represent the limiting 
envelope of trajectories, whose impingement points on the cowl surface 
now represent the limiting impact points. Any trajectories which start 
at (x o ,r o min<r o <r o max,0 Q ) will impinge on the surface region enclosed by 
the boundary curves defined by the inner and outer limiting impact 
points. 

4.6 Calculation of the Local Droplet Impingement Efficiency (8) 

One can now run a number of trajectories starting at an array of points 
within the region, r(r o min(0 Q ), r Q max(0 o )), in an orderly fashion to 
obtain the array of impingement points on the cowl surface (Figure 9). 
Let x (I.J) and 7 (n (I, J) denote the array of starting points from 
T(r min(0 ), r max(0 Q )) and the corresponding impingement point array 
on the surface, respectively. 
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impingement point array [x ( I , J ) ] 



Figure 9 - Illustration of the Starting Point Array [x Q (I,J)] and 
Impingement Point Array [x (I,J)] 


We can define the associated area elements; 


a Q (i,j) = area of the element formed by the four corner points, 
lx 0 (I,J). x 0 (I,J+l), X Q (I+l,J), X Q (I+1,J+1) } 

at the freestream constant-x plane, 

a m (i,j) = area of the element formed by the four corner points, 

u m (i,j), x m (i,j+i), x m (i+i,j), x m (r+i,j+i)} 

at the surface, 


where 

I = 

1 2 

1 

IMAX 

(IMAX 

= nmax) 


J = 

1 2 

1 • • • i 

JMAX 

(JMAX 

= nmax/2) 


i = 

1 2 

a » • • • i 

I MAX-1 




J = 

1 2 

JMAX- 1 




The local impingement efficiency (0), as in equation (1.0), can be 
approximated by 

A 

B(I c (i,j)) = n -.' U o (i ’ J1 (4-13) 

where x (i,j) = centroid location of a^ ( i , j ) . 

The unit direction vector, n , of V a for the general engine inlet 

attitude having a set of pitch (a), roll ($) and yaw (iji) angles with 

respect to the space coordinate axes can be obtained from the Euler 
rotation matrix applied to the unit vector ootained when the engine 
inlet body axes and the space axes coincide: 
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l"„(» n„(2) "J3)1 T 


‘ cosacosip 

sinotsin4>cosxp-cos4>sin4> 

si nacos4>cosiJ;+s i n<j>s i m|T 


cosasinij) 

s i n a s i n 4> s i n ij) + c o s 4> c o s 

sinacos<t>sin4i-sin<|>cos4> 


-sina 

cosasin^ 

cosacosefs 



cosacosil) 

cosasimj) 

-sina 


(4-14) 


Substituting (4-14) into (4-13), 

B(x c (i,j)) = cosacosil)a Q (i,j)/a m (i,j) (4-15) 

Equation (4-15) means that for a general orientation of the body with 
respect to the space coordinate system, we can choose the constant-x 
plane at the freestream as the plane of trajectory starting positions to 
compute the flux ratios. The projection of the freestream flux along 
the direction of 7 is accounted for by the factor (cosacosilO involving 
the pitch and yaw angles only. 

Thus the grid of 8 values can be computed numerically, using equation 
(4-15), at the centroids of the impingement point grid, x m (I,J). Unless 
the grids are dense, i.e., large IMAX and JMAX, the B distribution on 
the surface defined at x c (i,j)'s will not be smooth. Also, 7. ( i , j ) ' s 

are not particularly useful in organizing and presenting the 

computed B distribution on the surface because of its point 
function definition of B on 3-0 surfaces. 
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In order to obtain a more accurate 8 distribution and to express it 
as a function of the surface arc length(s) along a set of constant-9 
cowl contours, the following had to be done: 

(1) Select a desired 0 (constant-9 cowl contour along which 8 is 
to be determined) 

(2) From the coarse impingement centroid grid, 7. ( i , j ) , and the 
corresponding x o (I,J) grid computed, find the 1^ and (9-ray 
indices) that result in the condition 

9 C (I 1 .J) < 9 < 0 C (I 2 ,J) , J = 1,2, — , JMAX 

Determine KMAX: [9 (I ? ) - 0 (I^l/KMAX = 1° 

Find tangent trajectories for the rays, 

0 (K) = 0(Ij) + K ; K = 1,2 KMAX-1 

Compute a Q (k,j), a m (k,j), 8(x c (k,j)); k = 1,2 KMAX-2 

j = 1,2 JMAX-1 

Steps (1) thru (5) are repeated for all other 9's desired. 

Thus a much finer grid definition of 8 is obtained that encloses a 
particular constant-0 cowl contour desired. By parametric cubic 
splining of the contour curve, the arc length (s) of a point on the 
curve can be computed with respect to the cowl hi-lite position: 

(s = 0 is the hi-lite; s+ = outside; s- = inside). 

Moving the points along the cowl contour with selected As, 8 at a point 
along the contour can be interpolated from the surface grid cell corner 
values {8(x c (k,j)), B(x c (k , j-1) ) , 8(x c (k+l, j)) , 8(x c (k+l, j+1) ) } of the 
centroid grid that encloses the point in question (Figure 10). 


(3) 

(4) 

(5) 
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impingement point grid, x m (I,J) 


*”$•** impingement centroid grid, x c ( i , j ) 



constant-0 cowl contour 


x c (k,j+l ) 



(k+1 ,j+l ) 


B(Q) “= B(x c (k+1 ,j+l))-u-v + B(x c (k,j+l))-u*(l-v) + B(x c (k,j))-(l-u)«(l-v) 
+ 6(x c (k+l , j ) ) • ( 1 - u ) • v 

Figure 10 - Illustration of the Centroid Impingement Grid [x c (i,j)] 
and 0 Interpolation 
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4.7 Method for Solid Particles 


Most of the assumptions involved in the trajectory model are valid also 
for analyzing solid particles such as sand grains. One exception is the 
assumption of spherical shape - sand particles come in all shapes 
(irregular). If sand particles are characterized in terms of mass, 
however, one can still consider the shape to be spherical with an 
effective diameter characterized by their settling behavior in still air 
such as the terminal velocity and range. Stokes' diameter, for example, 
is the effective diameter of an equivalent spherical particle having the 
same mass and terminal velocity, based on the Stokes' drag, as the 
actual particle. Available data on sand particles of varying shapes 
( 43 , 44 | indicate that, as long as the size is not too large (<lmm), the 
settling behavior is about the same for particles having similar masses. 

The other consideration is the impact behavior - solid particles will 
bounce off the surface of impact. Shape of the particle will affect the 
bounce kinematics in a complicated way. Detailed analysis of such 
behavior is not worth pursuing, except to say that it will involve 
analyses of statistical nature. Some experimental data (451 are 
available that characterize the average behavior of solid particle 
kinematics at metallic surfaces of impact. 


The controlling kinematic parameter is the resiliency coefficients, r n 
and r t , relating the normal (U^) and tangential (U^) particle 
velocities after the impact to those before impact, and U^: 


U t2 " r t U tl 

The normal and tangential components of are 

U„1 = "(<Y"> 

U tl = n x (0^ x n) 


( 4 - 16 ) 


( 4 - 17 ) 
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where n = unit normal vector at the impact point (pointing outward 
from the surface) 

Combining (4-16) and (4-17), 


*n2 " 

-r n n( u i-n) 


U t2 ■ 

r^.n x (U^ x n) 


*2 = 

U n2 + °t2 = r t" X ^1 X ^ " r n^l^ 

(4-18) 


With known values of the resiliency coefficients, the particle velocity 
immediately after impact (U^) can be determined from equation (4-18). In 
the numerical integration of the trajectory equation, incorporation of the 
particle bounce mechanism amounts to restarting the integration process 
with the renewed initial condition, (x^, U^), at the point of impact. 


The impact point position (x m ) and the unit normal (n) are computed as 
described in the trajectory-bi-cubic patch intersection algorithm 
(Appendix A). 


Some available experimental data (45| indicate that the resiliency 
coefficients are functions of the impact incident angle, y , as well as 
the incident particle velocity magnitude, I U ^ I ; 


r n< IU l'- V 


't = r t (IU l'' V 


More research work in the measurement of these parameters are required in 
order to obtain an adequate empirical kinematic model for sand particles. 
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5.0 RESULTS AND DISCUSSIONS 


The basic numerical integration scheme of solving 3-D particle 
trajectories in potential flow of air was first checked by analyzing the 
water droplet impingement problem on a spherical body. This task was 
performed for a number of reasons: 

(1) Availability of the well documented NACA wind tunnel impingement 
data [Il|. 

(2) Well known incompressible potential flow solution in analytical 
form suitable for speedy code implementation. 

(3) Trajectory-surface intersection is easy to compute on the body. 

(4) 6 computation on the body is simple due to its axisymmetry. 

Having gained confidence from the axisymmetric analysis, the numerical 
scheme for a full-fledged 3-D impingement analysis was worked out, 
incorporating much of the code developed for the axisymmetric problem. 

5.1 ANALYSIS OF DROPLET IMPINGEMENT ON A SPHERE 

The potential flow velocity, V, can be expressed in terms of the 
Cartesian coordinates, (x,y,z): 

V(x,y,z) = [1 + ^r" 3 (L-3x 2 r -2 ) , -|xyr" 5 , |xzr" 5 | (5-1) 

2 2 2 2 
where r = x +y +z , 

(x, y,z) = (X,Y,Z)/R = non-dimensional field point 

R = radius of sphere located at (0,0,0) 
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The intersection between a trajectory segment, l(w), and the surface of 
2 2 2 

sphere, x +y +z = 1, is obtained by solving 

3 ? 

I (P n - + ^Lw r = 1 (5-2) 

i=l 

which is quadratic in w (w,L,|p and X are defined as in equation (A-34), 
Appendix A) . 


1 (w) = p + £Lw 


=> 


111 
X + y +2 


Calculation of the limiting envelope of trajectories follows the same 
procedure described in Section 4.5, except that the scan for 9 rays is 
not required due to the axisymmetry of the problem. However, several 
trajectories having different 9 starting values were checked to verify 
the axisymmetry in computed trajectories. 

The local impingement efficiency (8) takes the following simple analytic 
form (Figure 11): 

8(o>) = r Q dr Q /(rds) = - d ( r Q 2 /2 ! /d { cosuo | (5-3) 

2 2 2 

where r = x + y + z ; freestream starting radial position 
o ooo r 

of an impinging trajectory (dimensionless with R), 

a) = s/R; angle subtended by the impingement point at the 
origin (center of sphere). 
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Using the definition, oj = s/R, equation (5-3) is expressed in terms of 
the dimensionless arc length, s/R: 


B ( s/R ) = -d(r|/2]/d(cos(s/R) ] (5-4) 

For a number of impinging droplet trajectories (i=l,2,. . . ,n) equation 
(5-4) can be evaluated by a cubic spline of { (r^/2) ^ } as a function 
of (cos(s/R) . } . The cubic spline derivatives thus obtained determine 
B(s/R). 


5.1.1 Langmuir-0 Tunnel Droplet Size Distribution 

The numerical procedure was carried out utilizing the accepted tunnel 
droplet size distribution due to Langmuir (171. This distribution is a 
discretized plot of the cumulative LWC fraction versus the seven droplet 
sizes normalized to the mean volumetric diameter (MVD), as shown in 
Figure 12. 

For a tunnel cloud condition of a particular MVD, calculation of B 
involves weighting according to the multi-droplet size distribution: 

7 

B(s/R) = l c. B i (s/R) (5-5) 

i = l 

where B^ is the local impingement efficiency due to the droplets of 
diameter group d., c. is the fraction of the total LWC contributed by 
droplets in the diameter group d.. 

5.1.2 Results 

The following tunnel condition was used in the computation in order to 
compare with the NACA test results: 
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d/MVD 


FIGURE 12 - LANGMUIR - 0 MULTI-DROPLET SIZE DISTRIBUTION 
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R (sphere radius) 


9 inches 


V 

= 157 kts 

CO 

p 

s® 

= 28 in.Hg 

T 

= 50°F 

500 

MVO 

= 11.5 and 14.7 pm 


A plot of a typical water droplet trajectory is shown in Figure 13 for a 
14.7 urn diameter droplet. Also shown are the two potential flow 
streamlines (freestream starting position of the droplet trajectory was 
identical to that of the upper streamline). The droplet trajectory 
exhibits the inertial behavior through its departure from the high 
curvature portions of the upper streamline. 

Figure 14 shows a comparison of computed and test data of B vs. s/R for 
the two MVD tunnel clouds. Computed results are in good agreement with 
the experimental data, well within the reported experimental errors of 
10% in LWC and 6 % in MVO for the tunnel clouds measured. 

The significance of the Cunningham drag correction is indicated by the 
closer agreement, near the flow stagnation region, between the computed 
and test B values. This observed trend is understood in terms of the 
increased droplet impingement ■ by the smaller droplet size population of 
the Langmuir-0 droplet spectra, resulting from the appreciable 
Cunningham drag reduction for these droplets. The impinging droplet 

flux due to smaller droplets is more localized near the stagnation 
region than in the case of larger droplets because of the differences in 
their inertia. 

5.2 ANALYSIS OF DROPLET IMPINGEMENT ON A 3-0 ENGINE INLET 

Impingement analysis on the 737-300 prototype inlet was performed for 
the following flight conditions: 
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Figure 13 - Droplet Trajectory and Streamlines; Sphere 
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Figure 14 - Local Impingement Efficiency; Sphere 
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a = 0° and 25° 

V = 175 kts 

00 

T = 50°F 

Soo 

P* = 28.2 in. Hg 

M^.p = .625 (compressor face Mach no.) 

d = 30!im 

The a = 25° condition, close to the inlet design separation envelope, 
was selected as the 'worst* test case to thoroughly check the trajectory 
analysis code. These conditions represent a realistic full-power- 
takeoff flight situation, but are not representative of a typical icing 
condition. However, as to be seen from the computed results, these 

represent a severe flowfield situation in terms of the droplet 
trajectory computation, because of the high engine suction flow as well 
as the extreme angle of attack involved. 

The computational boundaries for the 3-0 full potential as well as the 
3-D trajectory code were 

|X ob " X CF I = 10r CF 

r max ” ^ r CF 

where X^ = freestream x-boundary of cylindrical computational 

doma i n 

X^p = compressor face x-boundary 
r^-p = inlet fan radius 

r = radial boundary for cylindrical computational 

max domain. 
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The computing time for each angle of attack case impingement analysis 

was about 5 min. on the CRAY-1S. 

5.2.1 Limiting Impingement Points 

The tangent impact points (*) are shown in Figures 15 and 16. The 

extent of droplet impingement is indicated by the portion of cowl 

surface bounded by these points. The 3-dimensional character of the 
flowfield and geometry causes the azimuthal variation of the tangent 

impact positions. It should be noted that no such variation in 6 is 
possible for an axisymmetric inlet at zero angle of attack, although 
9-variation of a different kind will result at non-zero angles of 
attack for axisymmetric inlets (non-axisymmetric flowfield). 

The following features are noted from Figures 15 and 16. 

(1) Wider extent of droplet impingement near 9 = 135° ('squashed' 
region) for both a = 0° and a = 25° cases. 

(2) For a = 0° case, all of the inner tangent impact points lie on the 

inner cowl surface; for a = 25°, a switch-over occurs at 9 = 110° 

beyond which the inner tangent impact points lie completely on the 
outer cowl surface. 

These observations can be understood in terms of (1) the increased 
exposure area to droplet impingement near the 'squashed' region due to 
the thickening of cowl cross-section near 9 = 135° and (2) the effect of 
high angle of attack causing the droplet impingement to occur more on 
the outer cowl surface for 9 > 90°. 

5.2.2 Droplet Trajectories on the Inlet Symmetry Plane 

The impinging droplet trajectories on the inlet symmetry plane are shown 
in Figures 17 thru 22. 
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The details of these near the upper (0 = 0°) and lower (0 = 180°) cowls 
are depicted in Figures 18 and 19 (a = 0°) and in Figures 21 and 22 
(a = 25°). 

The main feature noted in these figures is that the trajectories exhibit 
strong inward curvature as they approach the lip region. This is caused 
by the strong inlet suction flow typical of the ful 1 -power-takeoff 
setting. The combined effect of the high suction and high angle of 
attack is even more pronounced in Figures 20 and 22, where some extreme 
trajectory turn-arounds are seen near the lower cowl lip. 

5.2.3 Computed Local Impingement Efficiency Distributions 

Plots of computed B as a function of cowl contour arc length (s) 
at 0 = 0°, 45°, 90°, 135° and 180° are shown in Figure 23 for both 
angles of attack. Arc length (s) is the surface contour distance 
measured from the origin at the hi-lite, along a particular constant-0 
cowl contour curve. 

The following features are noted from the computed B curves; 

(1) 6 peaks broaden and decrease monotically from 9 = 0° to 0 = 135° 
for both angles of attack. 

(2) Zero angle of attack cases exhibit more rapid rise to the peak 8 
values compared to a = 25° cases, with the exception of curves for 
9 = 180°. 

(3) Double B peaks are observed for 9 = 180° (a = 25°), with a sharp 
maximum of 8 = .84 and a weak secondary peak of B = .25. 

Observation (1) can be explained in terms of the geometric character of 

the inlet; the progressive thickening of the cowl cross-sections 

from 9 = 0° to 0 = 135°, thus resulting in successively lower local 

droplet flux while increasing the extent of exposure to droplet 

impingement. 
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The second feature can be explained in that the distribution of 
impinging trajectories is more asymmetric about cowl hi-lites in a = 0° 
cases than a = 25°, as seen in Figures 18 and 21 for example. The 

strong inward curvature of trajectories adjacent to the inner tangent 
trajectory causes the location of the near normal impaction to be closer 
to the inner tangent impact point. For a = 25° cases the effect of high 
angle of attack partially offsets this trend, resulting in a more 
symmetric distribution. 

The presence of strong turn-arounds in trajectories (Figure 22) is 
responsible for the sharp primary 8 peak in observation (3), causing a 
localized region of high droplet flux. The weak secondary peak in 8 
occurs near the location of near normal impaction; local maxima in 8 is 
expected to occur whenever the impact velocity of droplet is nearly 
aligned with the surface normal vector, as shown in equation (1.0). 
Similar observations were also reported in 2-D cases, involving an 
airfoil with a leading edge ram air scoop [ 22 ) and ice-accreted airfoils 
[46|, where the trajectories undergo strong turn-arounds due to the 
abrupt and strong flow gradients present near the leading edge. 

The tails of the 8 curves are not plotted because these will require 
extrapolation of 8 in the present numerical scheme (see Figure 10); 
although the 8 values must go to zero at the tangent impact points (zero 

droplet flux since U-n = o there), the edge of the centroid 

grid { x” c ( i , j ) } is reached before getting to the edge of the impingement 
point grid {7 (n (I,J)}. 

Experimental impingement data for the 3-0 inlet analyzed are not 
presently available for comparison with the computed results. However, 
there is an on-going research program (Joint BMAC-Wichita State 
University) to obtain impingement data for the inlet analyzed as well as 
for several other geometries during the time period 1985-1986. Project 
director for this research program is Or. G.W. Zumwalt and the 

experimental 8 measurement involves a dye-water mixture spraying 

technique as well as laser reflectance spectroscopy. Tests will be 
conducted at the NASA-Lewis Icing Research Tunnel ( I RT ) under joint FAA 
and NASA sponsorship. 
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Figure 15 


- Limiting Impingement Points; a = 0° , 737-300 Prototype Inlet 
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igure 17 - Impinging Droplet Trajectories , a = 0 
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Figure 21 - Impinging Droplet Trajectories, Upper Cowl 
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Figure 23 - Computed Local Impingement Efficiency Distribution; 
737-300 Prototype Inlet, a = 0° and 25° 
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5.3 TRAJECTORY SIMULATION FOR A SOLID PARTICLE 


A sample trajectory simulation of a solid particle was performed using 
the flowfield about the 737-300 prototype inlet at 25° angle of attack. 
A simplistic bounce kinematic model assuming 50£ momentum loss along the 
normal and no loss along the tangential direction (r n = 0.5, r^ = 1.0 in 
equation (4-18)) was used to compute the particle velocity immediately 
after the impact. 

The particle was injected into the flowfield at about two fan radii 

ahead of the inlet with the initial velocity (cylindrical) of 

(U ,U ,U Q 1 = [0, -2 V , V 1. The ricochet trajectory of the particle is 
x r y c» <x> 

depicted in Figure 24. 

The sand separator efficiency of an engine inlet can be determined by 
tracing many such trajectories to compute the normalized freestream flux 
(or the flux through the surface of initial trajectory positions) of 
particles that correspond to the particle flux at the sand separator 
(scavenge) channel of the engine inlet. 
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FIGURE 24 - SAND PARTICLE TRAJECTORY TRACE 
737-300 PROTOTYPE INLET, a = 


6.0 CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER WORK 


Based on the results, findings and the physical assumptions of the 

present analysis method, the following conclusions can be made: 

(1) This investigation represents the first attempt in the analytical 
determination of the detailed water droplet impingement efficiency 
distribution on a 3-0 engine inlet surface. 

(2) Accurate surface definition of the 3-D surface is essential in the 
3-0 droplet trajectory/impingement analyses. 

(3) For the grid (mesh) definition of the flowfield, the flowfield 

accuracy as well as an accurate means of flow velocity 
interpolation are essential in computing accurate droplet 
trajectories. 

(4) Integration error control in the numerical integration of the 

trajectory equation is an important consideration in that it 
governs the computational efficiency as well as the accuracy of the 
computed trajectories with respect to the flowfield accuracy. 

(5) This analysis tool can easily be extended for problems involving 

solid particles, such as sand ingestion analyses. 

(6) Although experimental data is not yet available to directly verify 
the analysis results for the engine inlet analyzed, gQod agreement 
is obtained between the published test and computed results for an 
axisymmetric problem investigated. 

(7) Present analysis tool will not be appropriate for problems 

involving large concentration of water droplets (LWC > 10 gm/m ) 

or large droplets (d > 1000 yip). 
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Some recommendations for further work are: 


(1) Particle heat transfer equations should be incorporated in order to 
extend the present method to transonic/supersonic flows or flows in 
the engine compressor and turbine stages. 

(2) Reliable experimental drag data for large, non-spherical droplets 
is needed to extend the droplet size range of the method. 

(3) Effect of thick boundary layers on droplet trajectory should be 
studied for internal flow applications. 

(4) Improved trajectory computation can be achieved by using a body- 
fitted mesh definition of the flowfield and solving the trajectory 
equation in the transformed Cartesian mesh obtained through the 
metric of transformation for t^e particular body-fitted mesh 
chosen. 
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APPENDIX A 


Parametric Description of Curves and Surfaces 


A. 1.0 Parametric Curve Description 


A 3-dimensional curve segment can be represented parametrically by a 
polynomial expansion of the parameter variable u 
nmax 

r(u) = l a u n , u e [0,1] (A-l) 

n= 0 


n = 0,1,2,..., nmax 


where r = position vector along the curve, 

| A. L. 

a = n order vector coefficient, 
n 

The vectors corresponding to the parameter values of u=0 and u=l, i.e., 
r(0) and r(l), are the end point positions of the curve segment. 

The maximum power (nmax) of the parametric variable retained in the 
expansion determines the order of parametrization. Thus nmax=3 (or 
nmax=5) represents a cubic (or quintic) parametrization. The cubic and 
quintic representations are currently two of the most commonly used 
forms of analytic curve/surface construction techniques [39]. 


It will be shown later that these curve parametrization equations lead 
to surface parametric equations. 
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A. 1.1 Parametric Cubic Curve Representation 

From equation (A-l), the parametric cubic curve is 

— — — — 2 — 3 

r(u) = a Q + a^ u + a^ u + a^ u 

Equation (A-2) can be written in matrix form, 

7(u) = [ U ! [ A | T 

where [ U ] = [ 1 u u 2 u^ | 

[ A | = ( a Q a 2 a 3 | 

Differentiating equation (A-3) with respect to u, 

7 y (u) = [ 0 1 2 u 3u 2 | [ A ] T 

The end point (u = 0, u = 1) quantities of equations (A-2) and (A 

7(0) =(10001 [A| T 

7(1) =[11111 [ A 1 T 

7 U (0) = [01001 (A1 T 

7 U (1) =[0123| [A| T 

These can be expressed in matrix form, 

[7(0) 7(1) 7 U (0) 7 u ( 1 ) | T = [ N | [A 1 T 

where 1000 

[ N j = 1111 

0 10 0 
0 12 3 


(A-2) 


(A-3) 

(A— 4) 
-3) are 


(A-5) 
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The L.H.S. of equation (A-5) is just the vector array of the four curve 
end point conditions (positions and tangents) required to determine the 
four vector coefficients of the cubic parametrization. To solve for 
(A), equation (A-5) is to be inverted: 

( A j T = [ C | (7(0) 7(1) 7 U (0) 7 u (l) ] T (A-6) 


where 


(Cl = [NT 1 


0 

-3 

2 


0 0 0 
0 1 0 
3 -2 -1 

-2 1 r 


Substitution of (A-6) into (A-3) completes the cubic curve 
parametrization. 


r(u) = [U 1 [Cl [ r (0 ) r(l) r u (0) r u (l) J T 


(A- 7 ) 


Equation (A-7) implies that the cubic parametrization requires only the 
end point positions and slopes as inputs. 


Another approach of arriving at the same result is to express r(u) in 
terms of a set of four cubic polynomial blending functions, [ f - | ; 


r(u) = 

7(0) f L 

(u) + r u 

(0) f 2 

(u) + r( 1) f 

3 <u) * r u (L) 

f 4 (u) 

(A-8) 

= 

(^(u) 

f 2 (u) f 3 

(u) f 4 

(u)l (7(0) 

f u<°> 7(1) 7 U 

,(i)i T 


Equation 

(A-8) 

imposes 

the 

fol lowing 

constraints 

on the 

blending 

functions 

• 








f l 

i 

f i 

f 2 

f 2 f 3 

f 3 f 4 

f 4 ' 


u = 

0 1 

0 

0 

0 0 

1 0 

0 

(A— 9 ) 

u = 

1 0 

0 

L 

0 0 

0 0 

1 
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The cubic blending functions are 
4 

; m = 1, 2, 3, 4 


f_(u ) 
m 


y k 

~ L b u 
. mn 
n=l 


(A- 10 ) 


Coefficients b are determined easily from (A-9) and (A- 10 ) : 

mn 


m = 1; 


1 

0 

0 

0 


= b 


11 


b ll + b 12 + b 13 
b 12 


’14 


b 12 + 2b 13 + 3b 14 


'11 
’ 12 
*13 
* 14 


1 

0 

-3 

2 


Repeating the process for the remaining three, we have 


f x (u) = 1 

f 2 ( u > = 

f 3 (u) - 

f 4 (u) = 


-3u 2 + 2u 3 
3u 2 - 2u 3 
-2u 2 + 3u 3 
-u 2 + u 3 


(A-U) 


In matrix form, (A- 1 1 ) becomes 


[ fj(u) f 2 (u ) f 3 (u) f 4 (u) I = [ 1 u u‘ 


u 3 I 


1 

0 

-3 

2 


0 

0 

3 

-2 


0 

1 

-2 

1 


0 

0 

1 

1 


= [ u I ( c 


(A- 12 ) 


Substituting (A-12) into (A-8) , we see that the same parametric equation 
(A-7) results. 
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A. 2.0 Parametric Surface Description 


The parametric curve formulation can be extended to describe a surface 
element (patch) by allowing the vector coefficients in equation (A-l) to 
be functions of a second parameter. 

A. 2.1 Parametric Bi-cubic Surface Representation 

Expressing the vector coefficients {"a. ] of equation (A-2) as functions 
of a second parameter v, we have 

— — — — 2 — 3 

r(u,v) = a Q (v) + a^(v)u + a 2 (v)u + a 3 (v)u (A-13) 

Since equations (A-2) and(A-8) are equivalent, i.e., a\'s are simply the 
linear combinations of the curve end point position and tangent vectors, 
one can also introduce the parameter v into the equation (A-8) to 
describe a patch; 

r(u,v) = r(0,v)f 1 (u) + r u (0,v)f 2 (u) + r(l,v)f 3 (u) + ? u (l,v)f 4 (u) (A-14) 

If we had arrived at the curve parametric equation using the parameter 
variable name v instead of u and then introduced the second parameter u, 
we would have obtained 

r(u,v) = r(u,0)f 1 (v) + r y (u ,0) f ? (v) +r(u,l)f 3 (v) + ? v (u,l)f 4 (v) . (A- 15) 

Inspection of equations (A-14) and (A- 15) shows that each form uses a 
set of single parameter blending functions and hence is not symmetric in 
u and v. Physical meaning of these equations can be shown from Figure 
Al. 
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The parametric form of (A-14) is expressed in terms of the variable end 
points (7(0, v) and 7(1, v)) and end point tangents (7 u (0,v) and 7 (l,v)). 
Therefore it represents a u-parameter family of curves spanning the 
patch as the curve end conditions are continuously changed by sweeping v 
in the interval, [0,11 - Likewise, equation (A- 15) represents a v- 
parameter family of curves spanning the patch. In order to obtain a bi- 
cubic surface parametrization symmetric in u and v, the two equations, 

(A-14) and (A-15), must be combined. 

From equation (A-L4) 

7(u,0) = 7(0, 0)f L (u) + 7(l,0)f 2 (u) + 7 u (0,0)f 3 (u) + 7 u (l,0)f 4 (u) (A-16) 

7(u,l) = 7(0,l)f 1 (u) + 7(l,l)f 2 (u) + 7 u (0,l)f 3 (u) + 7 u (l,l)f 4 (u) (A- 17) 

Taking partials of (A-14) with respect to v 

r y (u,v) = 7 v (o,v)f 1 (u) + 7 y (l,v)f 2 (u) + 7 uv (0,v)f 3 (u) + 7 yv (l.v)f 4 (u) 

from which 

7 v (u,0) = 7 v (0,0)f 1 (u) + 7 v (l,0)f 2 (u) + 7 uv (0,0)f 3 (u) + 7 uv (l,0)f 4 (u) (A-18) 

7 v (u,l) = 7 y (0,l)f 1 (u) + 7 v (l,l)f 2 (u) + 7 uv (0,l)f 3 (u) + 7 uv (l,l)f 4 (u) (A- 19) 

Substituting equations (A-16) thru (A-19) into (A-15) and re-expressing each 
term on the R.H.S. of (A-15) in matrix form. 
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r(u,l)f 2 (v) = [ r(0,l)f 1 (u)+r(l,l)f 2 (u)+r u (0J)f 3 (u)+r u (l,l)f 4 (u)]f 2 (v) 


= [ ^(u) f 2 (u) f 3 (u) f 4 (u)] 

'0 F(0, 1 ) 0 o' 



0 F(1 ,1) 0 0 


f 2 ( y ) 


o F u (o,i) 0 0 


f 3 (v) 


_0 F u (l,l) 0 0_ 




(A-21) 


F v (u,0)f 3 (v)= [r v (0,0)f 1 (u)+r v (l,0)f 2 (u)+r uv (0,0)f 3 (u)+r uv (l,0)f 4 (u)]f 3 (v) 


[f^u) f 2 ( u ) f 3 ( u ) f 4 (u)] 

‘o o F v (o,o) o' 




o o F V ( 1 , 0 ) o 


f 2 (v) 


0 0 F uv (0,0) 0 


f 3 (v) 


_o o F UV (1 ,0) 0_ 




(A-22) 


F v (u,1)f 4 (v)= [r v (0,l)f 1 (u)+r v (l f l)f 2 (u)+r ijv (0,l)f 3 (u)+r uv (l.l)f 4 (u)]f 4 (v) 


= [^(u) f 2 (u) f 3 (u) f 4 (u)] 


0 0 0 F v (0,l)' 


’^(v) 

000 F v (l,l) 


f 2 (v) 

0 0 0 F UV (0,1) 


f 3 (*i 

_o o o F uv o,n 


[f 4 (v)J 


( A- 2 3 ) 


Adding equations (A-20) thru (A-23), equation (A- 15) becomes 


r (u ,v) = [f^u) f 2 (u) f 3 (u) f 4 (u)][Q][f 1 (v) f 2 ( v ) f 3 (v) f 4 (v)] T (A-24) 
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where the patch boundary matrix [ Q ] is given by 


[Q] s 


r(0,0) 

r(0,l) 

r v (0,0) 

F v<0.')' 

r(l,0) 

F(1 ,1) 

%('.») 

F v c.') 

r u (0,0) 

F(0,l) 

F uv<°'°> 

F uv (°.') 



F U v<>-°> 



(A-25) 


Using the relationship ( A- 12 ) the blending function arrays can be 
expressed as 

[f 1 (u)f 2 (u)f 3 (u)f 4 (u)| = fU) (C| 

|f 1 (v)f 2 (v)f 3 (.)f 4 (v)| T • (C| T (V| T 

where [U|, [V|, and [Cj are defined as before. 

Thus, the bi-cubic surface parametric equation is 

r(u,v) = [U] [Cl [Q | [C1 T [V] T (A-26) 

Looking at the array elements of [Q| the bi-cubic surface 

oarametrization requires the four vector quantities at each of the four 
corners of the patch element; 


74 



r = (position vector) 

r e (tangent vector along constant-v curve) 
r^ = (tangent vector along constant-u curve) 

? uv = (twist vector) 

Equation (A-26) is in a compact matrix form, well suited for numerical 
programming purposes. 

A. 2. 2 Geometrical Properties of Bi-cubic Surface Parametrization 

From the bi-cubic parametric patch equation (A-26) many geometrical 
quantities can be obtained in analytic forms. 

The unit normal vector (n) at a point r(u Q ,v o ) on the surface: 


n 


r (u„,v 1 x r (u .v ) 
u o o v v o o 

l r (u„ ) x F (u ,vj 
1 u o o v o o 


(A -27) 


(proper sign is to be chosen for particular application) 
Here, F u^ u o ,v o^ = CU' (u 0 )]rC][Q][C] T [V(v 0 )] T 

F v (u o’ v o ) = CU(u o )3[C][Q1CC] T [V(v o )] T 
[U'j 5 ^U] = [0 1 2u 3u l 2 ] 

[V'] = gfCV] = [0 1 2v 3v 2 ]. 
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The surface area (S) of a patch surface: 


1 i _ _ 

S = J / |r x r Idudv 

1 0 J 0 U V 1 

where (? u x F ),. = 

tfjjj = [U'][C][Qj][C] T [V] T 

(f£) k = [u][c][q k ][c] T [v'] T 


The volume (Vo 1 ) subtended by a patch surface at the origin: 


1 1 1 — — - 
Vo1 = 3 /„/„ r * (r u x r v )dudv 


where r*(r x r ) 
u v ' 


HI G i j k r i^ 3 u ^ j ^ 9 v ^ k 
i = 1 j = 1 k=l 


(A-28) 


(A-29) 
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A. 2. 3 Coordinate Transformations on Patch Surface 


Once a surface is parametrized, transformations, such as translation and 
rotation, are easily implemented without having to modify the parameter 
functions - only the vector quantities at the boundary corners of the 
patch need be transformed. If (R| denotes the particular transformation 
matrix and r* the position vector in the transformed system 


F*(u,v) = [U][C][Q*][C] T [V] T 


(A— 30 ) 


where components of the new patch boundary matrix [Q*] are 


3 

Q* = l R Q ; i = 1,2,3 
1 j=l 1J J 


( A- 31) 


Similarly derivative quantities are obtained using the same [Q*|; 


r u *(u,v) = [U'| [C| [Q*| [C| T [V| T ( A- 32 ) 
A. 3.0 Trajectory-Bicubic Patch Intersection 

To comoute the impingement point on the surface defined in terms of the 
bicubic parametric patch notation, equations (A-25) and (A-26), the 
geometric intersection between the trajectory line segment and the patch 
must be determined. 
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A. 3.1 Problem Definition 


At the required intersection point, the following must be satisfied: 


F.(u,v,w) = l.(w) - r . (u,v) * 0; i = 1,2,3 (A-33) 

Quantities T(w) and r(u,v) represent the position along the line 
(trajectory) segment and the position on the patch surface respectively: 


(line) 

where 


l.(w) = p. + C i wL (A-34) 

^ = (q” - pj/L = unit direction vector along the line, 

L e |p - q| = length of the line segment, 

w = normalized line length parameter ( 0^w<l ), 

(p,q) = line end position vectors. 


(patch) r.( u> ,) = [U]CC][Q i ][C3 T |:V] T 


A. 3.2 Numerical Method 

Solution of (A-33) is obtained by the Newton-Raphson technique for 
solving a system of non-linear algebraic equations: 

x'(n) = (x 1 , x 2 , x 3 ) s (u,v,w) 

F(n) = T[x 3 (n)] - rtx 1 (n),x 2 (n)] 

3F i 

Jij( p ) - tg x .^x(n) 

0 

3 

l J..(n)*A-(n) = - F.(n) ( A is solved by Gauss elimination ) 

j=l 1J J 1 

x(n+l ) = x(n ) + A(n ) 


( 1 ) 

( 2 ) 

(3) 

( 4 ) 
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Steps (1) thru (4) are iterated, each time updated with the correction 
vector A, until the equation (A-33) is satisfied within a specified 
tolerance. 

Since the solution is sought within the specified line segment and on 
the patch surface, the following constraints must be imposed on the 
independent variables: 

0 < x 1 (=u) < 1 
0 < x 2 (=v) < 1 
0 < x 3 (=w) < 1 

The constraint equations are satisfied by the transformation 

x. = s* /(I + s’) ; s i e ; i = 1,2,3 ( A -35) 

The Newton-Raphson steps can be modified accordingly in terms of the new 
variable s = (SpS^.s^): 


J?i( n ) = I 
1J k=l 


9F i 9x k , 

— - 6 .. 

9x k 3s j ° 


s(n ) 


■ J ij (n) ' 2 s j‘ ci + s j >' 2 


(A-36) 
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(A— 37) 


3 


l 0^,(n)-A^(n) = - F . [ s ( n ) ] 

j = l J 


s(n+l ) = s(n) + A*(n) 


(A- 38) 


Thus the transformation (A-35) modifies the Newton-Raphson Jacobian 
elements in a simple way, as shown in (A-36). The components of the 
unconstrained Jacobian elements are 

J 11 V w> • r , (“.*)] ■ - 3r ^ - U — = - [U'][C][0,][C] T [V] T 
J i2 = g“{ i,(«) - r ( (u,v)] = - 8r *‘ U ' V) • - CU3Cc:CQ i ]CC] T [V3 T 

J i3 V w) - = ?i L 

At the converged solution (u Q , v , w q ), the intersection point (7 m ) and 
the unit outward normal vector n(x m ) are computed by substitution of 
(u ,v Q ) into the patch parametric equations (A-26) and (A-27): 


H/ ( I H 2 )’' 2 
1 j = l J 

3 3 

j4k=i eiJk ’ Jj ’ (vv . ) ' Jk2<u °' v ° ) 


(x m>i * 

<">i ■ 

H. 

1 
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APPE NDIX 8 


Derivation of the Trajectory Equation (3-1): 


dU./dt = C 0 (R V )-R V - (V i -U i )/(24P) - (l-o)gC6 i2 /V* ( B ' l > 

A vector diagram of the forces acting on the particle is shown in Figure 
Bl. 


u 

y 



Figure Bl - Forces Acting on Particle and Velocities 
at Center of Particle 


Drag force, D, is given by 


°i = C Q (R v )-q-a-cos(i,n) 


(B-2) 


where R y 


q 

a 

V 

cos(i ,n) 


= p|V-U|d/y, 

= sphere drag coefficient, 

= 2 p | V-U | = velocity head experienced by the particle 

in the flowfield, 

2 

= Ttd /4 = projected area of spherical particle, 

= flow velocity, U = particle velocity, 

= (V-U)./|V-U| = direction cosine between the unit 

/S A 

vectors i ( // to i-axis) and n ( // to V-U). 
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Using the definitions, (B-2) becomes 


D i = 2‘pc D (R v ) - a -|v-u|-(v-u) i 


(8-3) 


Buoyancy and gravity forces are 


B i = pgx6 i2 


G.j = -p*gx6 ^ 2 ; x = nd /6 = volume of spherical particle 


(B— 4 ) 


where the gravity force is assumed to act along the -ve y (i=2) axis. 


Summing the forces ( ^ f = ma ) , 


dU./dt = ioC D (R v )-(a/x)*|V-UKV-U). - (l-o)g5 i2 


(8-5) 


2 

Introducing the' inertia parameter, P=p*d V ro /(18yC) , 

where C is the characteristic dimension of the boundary surface, the 
coefficient of j y - U | ( V - U ) - term in (8-5) can be 

written as 

\ aC D (R v )(a/x) = pC D ( R y ) - ( rrd 2 / 4 ) - ( 2f^rd 3 /6 ) _1 = 3pC D (R v ) • (4p*d) _1 

3pC d (R ) 

4 ( P*d Vqo \ ^ 1 8yC ^ 

18 yC V d 

C n (R ) p V d 
D v <» 

~24PC y 

W R 


24PC 
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where 


R 5 pVood/u ( Reynolds number based on d ). 


Equation (B-5) thus becomes 


du.j/dt = 


C 0 (R v ) -R* |V - U| (V - U). 


24 PC 


- g(l -o)6. 


i; 


Non-dimensional izing (B-6) 


d[ yu 

d[ t/(C /Vj] 


C n (R ) *R -(V./V -U. /V ) 

D v V' V v 1 oo 1 oo' 

24 P 


- (1 


which is 


dlK/dt 


C n (R )-R *(V.-U.) 
D v v ; v v i v 

24P 


- 0 - o)gC6 i2 /V^ 


where U- = U i /V cD , V. = V./V^ and t 2 t/(C/Vj. 


(B-6) 


-o)gC6 i2 /v2 
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APPENDIX C 


Table of Sphere Drag Coefficient (Cp) : 

# * ?< * 


R 

8TND.C 

C (STGHES) 

C (LANG) 

C (GAUVIN) 

C (PUTNAM) 

C (NORMENT 


D 

D 

D 

D 

D 

D 

0.05 

434.3200 

480. 0000 

494. 3258 

48*7. 1944 

490.3577 

436. 5822 

0. 10 

244.3200 

240.0000 

251.0861 

247.4012 

243.6177 

244.9004 

0.20 

124.4400 

120.0000 

128.5796 

125.9577 

126.3399 

124.0419 

0.40 

64. 3800 

60.0000 

66.6405 

64. 7953 

65.4238 

63.5772 

0.60 

44. 1200 

40 . 0000 

45.7168 

44 . 2242 

44.7425 

43.3909 

0.30 

34.2600 

30. 0000 

35. 1407 

33.8604 

34.3089 

33.2742 

1.00 

23.2240 

24.0000 

28. 7342 

27 . 6000 

23.0000 

27. 1854 

1.20 

24.0200 

20.0000 

24.4263 

23.4003 

23. 7641 

23.1104 

1.40 

21.0000 

17. 1429 

21.3245 

20. 3830 

20.7185 

20. 1862 

1.60 

18.7200 

15.0000 

13.9808 

13. 1075 

13.4200 

17.9314 

1.30 

16.8*733 

13. 3333 

17. 1450 

16.3234 

16.6216 

16.2560 

2.00 

15. 4200 

12. 0000 

15.6666 

14.3979 

15. 1748 

14.3662 

2.50 

12.7872 

9.6000 

12.9774 

12.3024 

12.5472 

12.3317 

3.00 

10.9*720 

3.0000 

11. 1533 

10.5525 

10.7734 

10.6027 

-3.50 

*7.6823 

6.8571 

9. 8414 

9.2894 

9.4917 

9.5507 

4.00 

3.6320 

6.0000 

3.3414 

3. 3327 

3.5193 

3.7009 

5.00 

7.2624 

4.8000 

7.4180 

6.9753 

7. 1392 

7.3395 

6.00 

6.2830 

4.0000 

6.4433 

6.0547 

6.2013 

6.4363 

8.00 

5.0340 

3.0000 

5.2042 

4.3777 

5.0000 

• 5. 1559 

10.00 

4.2763 

2.4000 

4.4313 

4. 1511 

4.2566 

4.3399 

12.00 

3.3020 

2.0000 

3.9013 

3.6539 

3.7472 

3.7759 

14.00 

3.4440 

1.7143 

3.5121 

3.2903 

3. 3739 

3.3631 

16.00 

3. 1635 

1 . 5000 

3.2128 

3.0115 

3.0874 

3.0479 

13.00 

2.9307 

1 . 3333 

2.9747 

2.7901 

2.8596 

2.7994 

20.00 

2.7492 

1 . 2000 

2.7801 

2.6095 

2.6736 

2.5934 

25.00 

2. 3*394 

0.9600 

2.4181 

2.2745 

2.3280 

2.2314 

30.00 

2. 1384 

0.8000 

2. 1659 

2.0415 

2.0873 

1 . 9327 

35.00 

1 . 9550 

0.6357 

1. 9785 

1 . 8683 

1 . 9086 

1 . 8027 

40.00 

1 . S078 

0 . 6000 

1 . 332*7 

1 . 7346 

1 . 7696 

1.6661 

50 . 00 

1.5970 

0.4800 

1.6195 

1.5381 

1 . 5653 

1 . 4720 

60 . 00 

1 . 4400 

0.4000 

1 . 4689 

1 . 3994 

1.4217 

1 . 3400 

30. 00 

1.2330 

0. 3000 

1. 2674 

1.2133 

1 . 2283 

1. 1699 

100. 00 

1 . 1016 

0. 2400 

1. 1363 

1.0917 

1. 1018 

1 . 0628 

120.00 

1 . 0020 

0.2000 

1 . 0427 

1 . 0045 

1.0110 

0.9375 

1 40 . 00 

0 . 9257 

0. 1714 

0. 9719 

0.9380 

0.9413 

0.9304 

160.00 

0.3640 

0. 1500 

0.9160 

0.3852 

0.8863 

0. 8343 

180.00 

0.3213 

0. 1333 

0. 3704 

0.3419 

0.3413 

0.8467 

200.00 

0. 7324 

0. 1200 

0. 3325 

0.3056 

0.3040 

0.3140 

250.00 

0.7085 

0.0*760 

0,7593 

0.7354 

0.7310 

0.7473 

300.00 

0. 660*3 

0.0800 

0.7075 

0.6839 

0.6775 

0.6934 

350.00 

0.6171 

0.0636 

0. 6676 

0.6440 

0.6362 

0.6469 

400 . 00 

0.53*72 

0.0600 

0.6359 

0.6119 

0.6029 

0. 6094 

500 . 00 

0.5501 

0.0480 

0.5885 

0.5627 

0.5520 

0.5627 

600.00 

0.5133 

0.0400 

0.5543 

0.5261 

0.5143 

0.5315 

300 . 00 

0.4743 

0.0300 

0.5077 

0.4743 

0.4609 

0,4926 

1000. 00 

0. 4469 

0.0240 

0. 4771 

0.4383 

0.4240 

0.4692 
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APPENDIX D 


Comparison between the Computed 3-D Full Potential CFD Flow Data 
and Wind Tunnel Test Data ( from NASA CR- 351 4 and Boeing Document 
06-49848): 


Wing pressures at: 



X 


Pylon pressures at 



Figure D1 


NASA Wing-Pylon-Nacelle Test Model; Location of 
Pressure Measurements Used for Analysis Comparison 
(from NASA CR-3514) 


# .63 1.38 

C = 24(1 + . 197R + 2.6E-04R )/R ; REFERENCED 3 

D 

* . 687 

C = 24(1 + . 1 5R >/R ; REFERENCED 3 

D 

S< 2/3 

C = 24(1 + 1/6R )/R ; REFERENCEC48 3 

D 

*2 2 3 

C R = 24.167R + 3.254R - . 23564R , . 05 £ R £l 3.04726 
D 

2 3 

= -23.339 + 38.969R + . 73204R - 5.6084E-04R , 3.04726 < R < 377.566 

2 

= 93.462R + . 37576R , 377.566 < R ; REFERENCE C49 3 
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Surface Mach number 



Test data: □ Upper surface 

<5> Lower surface 


Figure D2 - Wing Surface Mach Number Distribution, Wing-Pylon-Nacel le 
Model; M^ = 0.2, a = 5°, D = 4.5" 

(from NASA CR-3514) 
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lUrface Mac 


Lower row 


Upper row 



x inches 


Test data Upper row 


0 Lower row 

Figure 03 - Pylon Surface Mach 
Nacelle Model ; M 

CO 

Number Distribution, Wing-Pylon 
= 0.2, a = 5°, D = 4.5" 

(from NASA CR-3514) 











1.05 

.00 • 

-II 

X inches 

Test data £> 90-0° 

O 180 0° 

£> 270 0° 

Figure 04 - Nacelle Surface Mach Number Distribution, Wing-Pylon 
Nacelle Model; = 0.2, a = 5° , 0 = 4.5" 

(from NASA CR-3514) 





Surface Mach number 



Test data CD Upper surface 
O Lower surface 


Figure 05 - Wing Surface Mach Number Distribution, Wing-Pylon- 

Nacelle Model; M^ = 0. 6 , a = 3° , D = 4.5" 

(from NASA CR-3514) 
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Surface Mach number 



X inches 

Test data o Upper row 
O Lower row 


Figure D6 - Pylon Surface Mach Number Distribution, Wing-Pylon- 
Nacelle Model; M ro = 0. 6 , a = 3° , D = 4. 5" 

(from NASA CR-3514) 
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Surface Mach number 


1 4 


1 2 


1 0 



x inches 

Test data © 90.0° 

© 180 0° 

A 270.0° 

Figure D7 - Nacelle Surface Mach Number Distribution, Wing-Pylon-Nacelle 
Model; M = 0. 6 , a = 3° , D = 4. 5" 

CO 

(from NASA CR-3514) 
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-O 

73 
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ENGINE <£ o 
INLET £ 

0.5 



Figure D8 - Computational Hesh in the Vicinity of a Commercial - 
Transport Turbofan-Engine Type Inlet 
(from NASA CR - 35 14) 
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U) 
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Figure D9 - Cowl Surface Mach Number Distribution for a Corronerci al Turbofan-Engine 
Type Inlet; a = 25° 

(from NASA CR-3514) 
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